
if (!requireNamespace("doFuture", quietly = TRUE)) install.packages("doFuture")
if (!requireNamespace("igraph", quietly = TRUE)) install.packages("igraph")
if (!requireNamespace("readxl", quietly = TRUE)) install.packages("readxl")
if (!requireNamespace("Matrix", quietly = TRUE)) install.packages("Matrix")
if (!requireNamespace("progressr", quietly = TRUE)) install.packages("progressr")
if (!requireNamespace("furrr", quietly = TRUE)) install.packages("furrr")

library(igraph)
library(readxl)
library(Matrix)
library(progressr)
library(future)
library(furrr)
library(doFuture)

# 启用全局进度条
handlers(global = TRUE)
handlers(handler_progress(format = "[:bar] :percent :message"))

### 1. 加载PPI网络数据（预处理） ####
cat("步骤1: 加载PPI网络数据...\n")
file_path <- "41467_2019_9186_MOESM3_ESM.xlsx"
ppi_data <- read_excel(file_path, sheet = 1)
ppi_network <- graph_from_data_frame(ppi_data[, c("Protein_A_Entrez_ID", "Protein_B_Entrez_ID")], 
                                     directed = FALSE)

# 提取最大连通分量
cat("提取最大连通分量...\n")
comp <- components(ppi_network)
largest_comp_index <- which.max(comp$csize)
ppi_network <- induced_subgraph(ppi_network, which(comp$membership == largest_comp_index))

### 2. 读取药物靶点和疾病基因文件 ####
cat("步骤2: 读取药物靶点和疾病基因文件...\n")
drug_targets <- read.csv("成分-靶点-ledebouriellol.csv")
disease_genes <- read.csv("疾病-靶点.csv")
unique_disease_genes_ids <- unique(disease_genes$entrez_id)

# 仅保留网络中的基因
valid_nodes <- V(ppi_network)$name
drug_targets <- drug_targets[drug_targets$entrez_id %in% valid_nodes, ]
unique_disease_genes_ids <- unique(unique_disease_genes_ids[unique_disease_genes_ids %in% valid_nodes])

cat("有效药物靶点:", nrow(drug_targets), "\n")
cat("有效疾病基因:", length(unique_disease_genes_ids), "\n")

### 3. 预计算距离矩阵和节点度 ####
cat("步骤3: 预计算疾病基因距离矩阵和节点度...\n")
disease_nodes <- unique_disease_genes_ids
all_nodes <- V(ppi_network)$name

# 预计算所有节点度
node_degrees <- degree(ppi_network)
names(node_degrees) <- all_nodes

# 计算距离矩阵（仅疾病基因）
distance_list <- vector("list", length = length(disease_nodes))

# 使用基础进度条
pb <- txtProgressBar(min = 0, max = length(disease_nodes), style = 3)
for (i in seq_along(disease_nodes)) {
  disease_node <- disease_nodes[i]
  v_index <- which(all_nodes == disease_node)
  
  dists <- distances(ppi_network, v = v_index, to = V(ppi_network), mode = "all")
  dists[is.infinite(dists)] <- NA
  
  distance_list[[i]] <- as.vector(dists)
  setTxtProgressBar(pb, i)
}
close(pb)

# 创建稀疏距离矩阵
sparse_dist <- do.call(cbind, distance_list)
rownames(sparse_dist) <- all_nodes
colnames(sparse_dist) <- disease_nodes
sparse_dist <- as(sparse_dist, "sparseMatrix")
rm(distance_list)
gc()

### 4. 优化距离计算函数 ###
calculate_ds_fast <- function(targets, sparse_dist) {
  valid_targets <- targets[targets %in% rownames(sparse_dist)]
  
  if (length(valid_targets) == 0) return(NA)
  
  # 获取目标节点在矩阵中的行索引
  row_indices <- match(valid_targets, rownames(sparse_dist))
  
  # 计算每个靶点的最小疾病距离
  min_dists <- apply(sparse_dist[row_indices, , drop = FALSE], 1, min, na.rm = TRUE)
  
  # 移除Inf值（当某行全为NA时min会返回Inf）
  valid_min_dists <- min_dists[is.finite(min_dists)]
  if (length(valid_min_dists) == 0) return(NA)
  
  mean(valid_min_dists)
}

### 5. 优化网络邻近度计算###
compute_network_proximity_fast <- function(targets, sparse_dist, node_degrees, n_random = 1000) {
  # 计算实际距离
  d_real <- calculate_ds_fast(targets, sparse_dist)
  if (is.na(d_real)) return(NA)
  
  all_nodes <- rownames(sparse_dist)
  valid_targets <- intersect(targets, all_nodes)
  
  # 获取有效靶点的度
  target_degrees <- node_degrees[valid_targets]
  
  # 创建度匹配的随机靶点池
  random_distances <- numeric(n_random)
  valid_count <- 0
  max_attempts <- n_random * 3  # 设置最大尝试次数
  
  # 度匹配抽样函数
  generate_degree_matched_sample <- function() {
    # 度匹配抽样：为每个靶点寻找度相近的替代节点
    sapply(target_degrees, function(deg) {
      # 寻找度在±15%范围内的节点
      deg_range <- c(max(0.85 * deg, 0), 1.15 * deg)
      candidates <- names(node_degrees)[node_degrees >= deg_range[1] & node_degrees <= deg_range[2]]
      
      # 如果候选太少，放宽范围
      if (length(candidates) < 10) {
        candidates <- names(node_degrees)[
          node_degrees >= max(1, 0.7 * deg) & node_degrees <= 1.3 * deg
        ]
      }
      
      if (length(candidates) > 0) sample(candidates, 1) else NA
    })
  }
  
  # 生成随机样本
  attempts <- 0
  while (valid_count < n_random && attempts < max_attempts) {
    attempts <- attempts + 1
    
    random_targets <- generate_degree_matched_sample()
    random_targets <- na.omit(random_targets)
    
    if (length(random_targets) > 0) {
      d_random <- calculate_ds_fast(random_targets, sparse_dist)
      if (!is.na(d_random)) {
        valid_count <- valid_count + 1
        random_distances[valid_count] <- d_random
      }
    }
  }
  
  # 如果有效样本不足，使用现有样本继续
  if (valid_count < n_random) {
    warning(paste("只生成了", valid_count, "个有效随机样本，而不是", n_random))
    random_distances <- random_distances[1:valid_count]
  }
  
  # 计算Z值
  mu_d <- mean(random_distances)
  sigma_d <- sd(random_distances)
  
  if (sigma_d < 1e-10 || is.na(sigma_d)) {
    return(ifelse(d_real < mu_d, -10, 10))
  }
  
  (d_real - mu_d) / sigma_d
}

### 6. 计算主流程 - 使用顺序计算避免并行问题 ###
drugs <- unique(drug_targets$chemical)
cat("步骤4: 计算药物网络邻近度...\n")
cat("药物数量:", length(drugs), "\n")

# 预计算所有药物靶点列表
drug_targets_list <- lapply(drugs, function(d) {
  unique(drug_targets$entrez_id[drug_targets$chemical == d])
})
names(drug_targets_list) <- drugs

# 开始计算
start_time <- Sys.time()
results <- data.frame()

# 使用顺序计算
for (i in seq_along(drugs)) {
  drug <- drugs[i]
  cat(sprintf("处理药物 %d/%d: %s\n", i, length(drugs), drug))
  
  targets <- drug_targets_list[[drug]]
  
  d_s <- calculate_ds_fast(targets, sparse_dist)
  Z_score <- compute_network_proximity_fast(targets, sparse_dist, node_degrees)
  P_value <- if(!is.na(Z_score)) 2 * pnorm(-abs(Z_score)) else NA
  
  results <- rbind(results, data.frame(drug, d_s, Z_score, P_value))
}

cat(sprintf("计算完成! 总耗时: %.1f秒\n", difftime(Sys.time(), start_time, units = "secs")))

### 7. 结果处理与输出 ###
cat("步骤5: 验证结果并输出...\n")
results$significance <- ifelse(results$Z_score < 0 & results$P_value < 0.05, 
                               "显著有效", "不显著")
results <- results[order(results$Z_score), ]

# 打印和保存结果
print(head(results, 10))
write.csv(results, "network_proximity_results.csv", row.names = FALSE)

# 可视化
if (nrow(results) > 1) {
  hist(results$Z_score, main = "药物网络邻近度Z值分布", 
       xlab = "Z值", ylab = "药物数量", col = "lightblue")
  abline(v = 0, col = "red", lwd = 2)
  text(0, max(hist(results$Z_score, plot = FALSE)$counts)/2, "Z=0", pos = 4, col = "red")
} else {
  cat("只有一种药物，跳过直方图绘制\n")
}

# 显示总耗时
end_time <- Sys.time()
cat(sprintf("总耗时: %.1f分钟\n", difftime(end_time, start_time, units = "mins")))